#ifndef _AL_SSI_H
#define _AL_SSI_H

#include <cmath>
#include <fstream>
#include <vector>
extern "C"{
#include <matrix2.h>	// header for mesch12b
}
#include "../al_structs.h"

/*static double LegendreTable[8][9]=
{
{0.0},
{0.0, 1.0},
{0.0},
};*/

template<typename ITER, int M, int N>
class TSphericalInterpolation
{
private:
	
	size_t inSize;
	VEC *vCi;
	ITER iterF;
	alArray2D *coefMI;
	alArray2D *LegendreTable;
	double dotP(const alSSIPoint &p1, const alSSIPoint &p2){
		return p1.x*p2.x+p1.y*p2.y+p1.z*p2.z;
	}
	double Gm(double x){
		double ret=0.0;
		for(size_t i=0;i<N;i++)
		{
			ret+=coefMI->getAt(M,i)*Legendre(i,x);
		}
		return ret;
	}
public:
	TSphericalInterpolation(const ITER first, const ITER last, alArray2D *coef, alArray2D *lt):iterF(first), coefMI(coef), LegendreTable(lt){
		MAT* matGMM=0;
		inSize=0;
		ITER start=first;
		while(start!=last)
			inSize++;
		matGMM = m_get(inSize, inSize);
		for(size_t i=0;i<inSize;i++){
			for(size_t j=0;j<inSize;j++)
			{
				if(j<i)
					matGMM->me[i][j] = matGMM->me[j][i];
				else
				{
					matGMM->me[i][j]=Gm(dotP(*(first+i),*(first+j)));
				}
			}
		}//for i
		PERM	*pivot;
		VEC		*b;
		b = v_get(inSize);
		for(size_t i=0;i<inSize;i++)
			b->ve[i] = (*(first+i)).value;
		pivot = px_get(matGMM->m);
		LUfactor(matGMM,pivot);
		/* set values of b here */
		vCi = LUsolve(matGMM,pivot,b,VNULL);
		M_FREE(matGMM);
	}//func
	
	~TSphericalInterpolation(){
		V_FREE(vCi);
	}
	
	double GetInterpolatedValue(double xpos, double ypos, double zpos){
		alSSIPoint pt(xpos, ypos, zpos);
		double ret=0;
		for(size_t i=0;i<inSize;i++){
			double cosine = dotP(pt, *(iterF+i));
			ret+= vCi->ve[i]*Gm(cosine);
		}//for i
		return ret;
	}
	static void ReadCoefFromFile(alArray2D **coef, alArray2D **lt){
		(*lt) = alArray2D::read("_al_legendre.dat");
		(*coef) = alArray2D::read("_al_coef.dat");

		(*coef)->Dump();
		(*lt)->Dump();
	}
	double Legendre(unsigned int n,double x){
		double ret=0.0;
		for(size_t i=0;i<n+1;i++){
			//ret += LegendreTable[n][i]*
			ret = LegendreTable->getAt(n,i) + ret*x;
		}
		return ret;
	}
};

typedef TSphericalInterpolation< std::vector< alSSIPoint >::const_iterator , 3, 5>  SphericalInterpolationV ;
#endif